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Abstract 

Mathematical models are increasingly being used to understand complex biochem- 
ical systems, but we rarely know the consequences of choosing one model to another. 
Using algebraic techniques we study systematically the effects of intermediate (tran- 
sient) complexes and provide a simple, yet rigorous mathematical classification of all 
models obtained from a core model by including intermediates. Main examples in- 
clude enzymatic and post-translational modifications systems, where intermediates 
often are considered insignificant and neglected in a model. In each class, a math- 
ematically simple canonical model characterizes crucial dynamical properties, such as 
mono/multistationarity and stability of steady states, of all models in the class. Im- 
portantly, our results provide guidelines to the modeler in choosing between models 
and in distinguishing their properties. 

Systems biology aims to understand complex systems and to build mathematical models 
that are useful for inference and prediction. However, model building is rarely straightfor- 
ward and we typically seek a compromise between the simple and the accurate, shaped by 
our current knowledge of the system. Transient, or intermediate, complexes in biochemical 
reaction pathways are often ignored in models or grouped into a single or few complexes, 
either for reasons of simplicity or conceptual clarification, or because of lack of knowledge. 
For example, models of the MAPK signaling cascade vary considerably in the details of 
intermediates [10^ [T5] and intermediates are often ignored in models of ubiquitous phospho- 
relays and two-component systems jH|13j. However, inclusion/exclusion of intermediates 
might affect the properties of a model and thus influence our reasoning and conclusions. 
An example is sequestration of intermediates that can cause ultrasensitive behavior [IJ [22] . 
It is therefore important to understand the role and consequences of model choice and 
model uncertainty in modeling biochemical systems. 

Algebraic techniques have the potential to provide qualitative insight about models 
in the lack of quantitative information about kinetic parameters and species abundances 
[Hi [121 Q3S EI] • Building on previous work [HI [20] , we propose a mathematical framework 
to study the dynamical properties of models that differ in how intermediates are included. 
The most fundamental and crucial dynamical features are the number and stability of 
steady states. Here we study the capacity of models to exhibit different steady-state 
features, assuming the kinetic parameters are unknown. 

An intermediate is defined as a species in a reaction network that is created and disso- 
ciated in isolation, that is, it is produced in at least one reaction, consumed in at least one 
reaction and it does not directly interact with any other species. A core model is the min- 
imal reaction mechanism to be modeled. Each reaction yi — > yj in the core model consists 
of two core complexes yi,Vj, which consist of core species. An extension model contains 
the core complexes together with some intermediates that are not part of the core model, 
and such that the core model can be obtained from the extension model by collapsing all 
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reaction paths yi — )■ Y\ — > . . . — > — > yj, where Yi are intermediates, into a single reaction 
y { yj (Fig. 1). 

One concretization of Fig. 1 is the simple enzymatic mechanism obtained by choosing 
yi = So + E, y 2 = S\ + E, and y 3 = S 2 + E, where E is an enzyme and Sj is a substrate 
with j = 0,1,2 phosphorylated sites. The substrate So can be doubly phosphorylated 
sequentially via Si or directly (processively) . Under mass-action kinetics, the dynamics of 
Fig. ID is described by a polynomial system of ordinary differential equations (ODEs): 

[So] = ~h [S ] [E] , [Y] = k\ [S ] [E] + k 3 [Si] [E] - k 2 [Y] - k 4 [Y] , 

[Si] = -k 3 [S!][E\ + k 2 [Y], [E] = -hiSoM - US^E] + k 2 [Y] + k A [Y], (1) 

[S2] = h[Y], 

where fc* are rate constants, [X] denotes the concentration of X, and [X] the instantaneous 
change in [X] . In addition there are two conservation laws, 

S£tal °= [So] + [Si] + [S 2 ] + [Y] , and < al D ^ f [E] + [Y] , (2) 

that is, quantities that are conserved over time and determined by the initial concentrations. 
Conservation laws confine the dynamics to an invariant space given by S^tai an d -E'totaii 
and the dynamical analysis must be restricted to this space. Also the core model in Fig. 1A 
has two conservation laws, 

^otai= f [5 ] + [5i] + [5 2 ], and Ef ot ^ [E]. (3) 

The similarity between ([2]) and Q holds generally (for proofs, see supplementary text): 

Fact 1: The conservation laws in the core model are in one-to-one correspondence with 
the conservation laws in any extension model. The correspondence is obtained by adding 
a suitable linear combination of the [Y] 7 s to each conservation law of the core model. □ 
We next state two facts that allow us to relate the dynamics near steady states of the 
core and extension models to each other. At steady state [X] = for all species X. We let 
[y] denote the product of the species concentrations in complex y, for example, [2S] = [S] 2 
and [^o+.E'] = [5o][-E7]. Different extension models contain different intermediates, resulting 
in different steady-state equations. However the intermediates can always be eliminated 
from the equations: 

Fact 2: The steady-state concentrations of the intermediates Y in an extension model are 
given as linear sums [Y] = Yly ^Y,y[y] of the concentrations of the core species, [y] appears 
in the expression, that is, fj,y,y 7^ 0, if and only if there is a reaction path y — > . . . — > Y 
involving exclusively intermediates [3, [8] . □ 
The constant fiy,y is either zero or positive and depends only on the rate constants in 
the extension model. In ([!]), the equation [Y] = gives 

[Y] = m 1 [So][E]+m 3 [S 1 ][E], (4) 

where mi = k ^_ kA are reciprocal Michaelis-Menten constants [3]. If Q is substituted into 
([!]), we obtain a new ODE system: 

[So] = -h[S ][E], [S\] = -hm^S^E] + k 2mi [S ][E], 

[E] = 0, [S 2 ] = k Ami [So] [E] + fc 4 m 3 [Si] [E] , (5) 
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Figure 1. Core and extension models. A core model and three possible extensions 
that incorporate intermediates in different ways. All intermediates (capital Y) are produced 
and subsequently consumed and cannot be part of another complex. A. A core model with 
three core complexes (yi) and three irreversible reactions. B. Three transient intermediates 
are involved in the reactions. C. The complex y\ is involved in a reversible "dead-end" 
reaction with one intermediate. D. The complex y\ is converted into Y, which splits into 
2/2 or j/3, respectively (the former reversibly). The core complexes might take any form, e.g. 
yi could be 5-1- E or 2 A. The core model is obtained by collapsing reaction paths. Observe 
that the directionality of the arrows needs to be preserved. For instance, an extension of 
the reaction y\ — > y 2 cannot be ?/i ^ 7 ^ 1/2, because it would imply that yi — > y\ is in 
the core model. By adding arbitrarily many intermediates (e.g. y\ ^ Y\ ^ . . . ^ Yjt) we 
can create arbitrarily many extension models with the same core. 



which is a mass-action system for the core model in Fig. 1A with k\ = ^^l, ^2 = ^4^13 , 
and &3 = fojmi (because k\ = k^vtii -\-k/^m\). We say that the rate constants fc* are realized 
by k*. In some (unrealistic) extension models, not all choices of rate constants in the core 
model are realizable (supplementary text). The relation between the ODEs in Fig. 1A and 
ID holds generally for any pair of core and extension models: 

Fact 3: After substituting the expressions [Y] = J2 y ^Y,y[v] into the ODEs of the extension 
model, we obtain a mass-action system for the core model. □ 

Consequently, core and extension models fulfill the same steady-state equations for the 
core species (Facts 2 and 3) but conservation laws impose different constraints (Fact 1). 
Specifically, by inserting Q into ^ we obtain 

Cai = [So] + [Si] + [S 2 ] + m 1 [S ][E] + msiS^E], 

<tal = [E] + mi [So] [E] + m 3 [Si] [E] . (6) 

The steady states of the extension model solve ([5| and ([6]), while they solve ([5j and ^ in 
the core model. The number of solutions can differ substantially (supplementary text). 

In view of these observations, we provide a classification of the extension models based 
on the core complexes that contribute to the conservation laws (if any). In (pi), these 



Simplifying Biochemical Models 



4 



Class complex yi 


Class complexes y\ , y2 


Class complexes y 2 , j/3 


3 steady 
states 


2/2 

Y / \ 
]{/ \ 


3 steady 
states 


y l / \ 

11/ \ 


1 steady 
state 


Y 2 ^V2 

/ \ F i 
/ \ll 


2/1 > 2/3 


2/1 > 2/3 


2/1 > 2/3 


2/2 2/2 

Y 2 ^ T\ 7 \ 

T »7i \ 3 \ 
2/1 2/3 2/1 Y" 2 — > 2/3 


2/2 2/2 

/ v. TI 
Y x Y 2 Y 
r \ /* n 

2/i — > y 3 — > 2/3 2/i 2/3 


2/2 2/2 

/ \ Tl 

/ yl y v 

/ \\ / X 

2/1 > 2/3 2/1 2/3 



Figure 2. Steady-state classes. We consider the core model of Fig. 1 and its steady- 
state classes. The canonical model is the extension model with a dead-end reaction added 
for each class complex (upper right corner). Each class (except the class of the core model) 
has an infinite number of members and a few of these are shown. For the number of 
steady states we consider the model with y\ = Sq + E, y 2 = Si + E, y% = S 2 + E and 
dephosphorylation reactions S 2 — > Si — >■ So (not shown in the figure). The number of 
steady states in each class refers the maximal number of steady states that a model in the 
class can have for some choice of rate constants and conserved amounts. See supplementary 
text for details. 

complexes are Sq + E and Si + E. Two extension models belong to the same steady- 
state class if they share the same core complexes in the conservation laws. The complexes 
characterizing a steady-state class are called the class complexes. 

Fact 2 provides a graphical characterization of the class complexes of a steady-state 
class (Fig. 2). The class of the core model has no complexes. The extension model in 
Fig. 1C belongs to the steady-state class with class complex yi because there is only one 
path from a core complex to an intermediate: yi — > Y. Similarly, the extension models 
in Figs. IB and ID have class complexes 2/1)2/2- We conclude that Figs. IB and ID are 
in the same class with similar equations for the conservation laws, while the models in 
Figs. 1A and 1C are in different classes and have different equations. In this case, Fig. 1A 
has always one steady state for any choice of conserved amounts and rate constants, while 
Figs. IB-ID can be multistationary (supplementary text). 

Since class complexes characterize the steady-state classes, there is a finite number of 
classes, 2^, with N the number of core complexes (N = 3 in Fig. 1). We say that a class 
is smaller than another class if the second contains the class complexes of the first. The 
steady-state class of the core model is smaller than any other class. A canonical model of 
a class can be constructed by only considering dead-end reactions, y ^ Y (Fig. 2). We 
now come to the main results concerning the dynamical properties of the classes. 

Fact 4: If the core model has non-degeneratd^ steady states for some rate constants 

3 A steady state is said to be non-degenerate if the Jacobian of the ODE system evaluated at the steady 
state is non-singular (supplementary text). 
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"All" extension models are multistationary 



Ql: Is the core model 
multistationary? 



None of the models 
are multistationary 



Q2: Is the maximal 
canonical model multi- 
stationary? 



"All" extension models 
that include the model are 
multistationary 



Q3: Is a canonical model 
multistationary? 



None of the models in- 
cluded in the model are 
multistationary 



Figure 3. Modeler's decision tree. A decision tree to detect multistationary steady- 
state classes. "All" means that the model exhibits multistationarity as long as the rate 
constants of the core model can be realized by the extension model. Q3 must be checked 
for different canonical models as necessary. 



and conserved amounts, then any extension model that realizes the rate constants has at 
least N corresponding non-degenerate steady states for some rate constants and conserved 
amounts. These can be chosen such that the correspondence preserves unstable steady 
states with at least one eigenvalue with non-zero real part and asymptotical stability for 
hyperbolic steady states. Reciprocally, if the canonical model of a steady-state class has 
a maximum of N steady states for any rate constants and conserved amounts, then all 
extension models in the class, or in any smaller class, have at most N steady states. □ 

In particular, if the largest canonical model (with a dead-end at all core complexes) 
cannot be multistationary, then no extension model, including the core model, can be 
multistationary. By screening the canonical models for the possibility of multistationarity, 
the modeler can obtain a clear idea about the effects of intermediates. In Fig. 2, the steady- 
state class given by 2/2)2/3 does not have multiple steady states, hence the same holds for 
the class given by 2/2 or 2/3 alone (Fact 4). Existence of multistationarity in Figs. 1B-D 
is due to the non-linearity introduced by [yi] in the conservation laws, irrespectively of 
whether [2/2] and [2/3] are included. If there are no conservation laws then any steady state 
in the core model corresponds precisely to a steady state in the extension model (assuming 
rate constants are realizable; Facts 2 and 3). 

Our approach provides a simple graphical procedure to classify the extension models 
into a finite set of classes with common dynamical features, thereby elucidating the con- 
sequences of choosing a specific model. Fig. 3 shows a decision diagram that guides the 
modeler through a number of possibilities. Checks can be performed using different com- 
putational methods [21 [6] or by manually solving the system (a task that simplifies due 
to the simple form of the canonical models). Fig. 4 shows a biological application of the 
decision diagram in Fig. 3. 

Our work develops from the perspective of the model and clarifies the effects of inter- 
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Figure 4. Multistationarity in two-component systems. Example of an applica- 
tion of the decision tree in Fig. 3. Four models of two-components systems are considered. 
All models are core in the sense that they are not extension models of any smaller models. 
SK=sensor kinase; RR=response regulator; Ph=phosphatase; T=phosphatase; *=phos- 
phorylated (activated) state. (A) Basic phosphorelay mechanism: SK autophosphorylates 
and transfers the phosphate group to RR; a phosphatase dephosphorylates RR*. (B) Same 
as (A), in addition SK is bifunctional and dephosphorylates RR*. (C) Same as (A), in 
addition SK and RR form a complex Y. (D) Same as (C) with the addition of a phos- 
phatase T for SK*. Systems (B)-(D) are core models of the models considered in [HIE]. 
The models analyzed in ITT] are extension models belonging to multistationary classes 
(last column of the table) and hence display multistationarity. The answers to Q1-Q3 have 
been obtained using the CRN toolbox [5]. 

mediates in biochemical modeling. Simplifications are always applied in model building 
but generally on a case to case basis, motivated by biological assumptions. One example 
is the Quasi-Steady-State Approximation (QSSA), where equations of the form [Y] = 0, 
together with some (but not all) conservation laws, are used to eliminated species [31 118]. 
This results in a hybrid model between our core and extension models. Model simplifi- 
cation and choice must be pursued with great care as crucial dynamical properties might 
change radically by the inclusion of intermediates. 

Our analysis is a first step in quantifying the relationship between simple and complex 
models of the same system, and in using simple models to predict properties of complex 
systems. Our results can guide the modeler through the critical issue of choosing a model 
and in learning about model properties. As such the results are useful for interpretation of 
experimental data and for designing synthetic systems. We envisage that our techniques 
can be extended to other models than those defined by intermediates and can provide 
further insight into the nature of biochemical and other modeling [7J 116] . 



Simplifying Biochemical Models 



7 



References 

[1] N. Bluthgen, F. J. Bruggeman, S. Legewie, H. Herzel, H. V. Westerhoff, and B. N. 
Kholodenko. Effects of sequestration on signal transduction cascades. FEBS J., 
273:895-906, 2006. 

[2] C. Conradi, D. Flockerzi, J. Raisch, and J. Stelling. Subnetwork analysis re- 
veals dynamic features of complex (bio)chemical networks. Proc. Nat. Acad. Sci., 
104(49):19175-80, 2007. 

[3] A. Cornish-Bowden. Fundamentals of Enzyme Kinetics. Portland Press, London, 
third edition, 2004. 

[4] A. Csikasz-Nagy, L. Cardelli, and OS. Soyer. Response dynamics of phosphorelays 
suggest their potential utility in cell signalling. J R S Interface, 8(57):480-8, 2011. 

[5] P. Ellison, M. Feinberg, H. Ji, and D. Knight. Chemical reaction network toolbox, 
version 2.2. http://www.chbmeng.ohio-state.edu/ feinberg/crntwin/, 2012. 

[6] E. Feliu and C. Wiuf. Preclusion of switch behavior in reaction networks with mass- 
action kinetics. Appl. Math. Comput, 219:1449-1467, 2012. 

[7] E. Feliu and C. Wiuf. Variable elimination in chemical reaction networks with mass- 
action kinetics. SI AM J. Appl. Math., 72, 2012. 

[8] E. Feliu and C. Wiuf. Variable elimination in post-translational modification reaction 
networks with mass-action kinetics. J. Math. Biol., In press:DOI: 10.1007/s00285- 
012-0510-4, 2012. 

[9] H. A. Harrington, K. L. Ho, T. Thorne, and M. P. H. Stumpf. Parameter-free model 
discrimination criterion based on steady-state coplanarity. Proc. Natl. Acad. Sci. 
U.S.A., 2012. 

[10] C. Y. Huang and J. E. Ferrell. Ultrasensitivity in the mitogen-activated protein kinase 
cascade. Proc. Natl. Acad. Sci. U.S.A., 93:10078-10083, 1996. 

[11] O.A. Igoshin, R. Alves, and M.A. Savageau. Hysteretic and graded responses in 
bacterial two-component signal transduction. Mol. Microbiol, 68:11961215, 2008. 

[12] R.L. Karp, M. Perez Millan, T. Dasgupta, A. Dickenstein, and J. Gunawardena. 
Complex-linear invariants of biochemical networks. Journal of Theoretical Biology, 
311:130-138, 2012. 

[13] JR. Kim and KH. Cho. The multi-step phosphorelay mechanism of unorthodox two- 
component systems in e. coli realizes ultrasensitivity to stimuli while maintaining 
robustness to noises. Computational biology and chemistry, 30(6):438-44, 2006. 

[14] E. L. King and C. Altman. A schematic method of deriving the rate laws for enzyme- 
catalyzed reactions. J. Phys. Chem., 60:1375-1378, 1956. 



Simplifying Biochemical Models 



8 



[15] N. I. Markevich, J. B. Hoek, and B. N. Kholodenko. Signaling switches and testabil- 
ity arising from multisite phosphorylation in protein kinase cascades. J. Cell Biol., 
164:353-359, 2004. 

[16] S. Rao, A. van der Schaft, K. van Eunen, B. M. Bakke, and B. Jayawardhana. Model- 
order reduction of biochemical reaction networks. arxiv:1212.2438, 2013. 

[17] B. Salvado, E. Vilaprinyo, H. Karathia, A. Sorribas, and R. Alves. Two component 
systems: Physiological effect of a third component. PLoS ONE, 7(2):e31095, 2012. 

[18] L. Segal and M. Slemrod. The quasi-steady-state assumption: A case study in per- 
turbation. SIAM Review, 31:446-477, 1989. 

[19] G. Shinar and M. Feinberg. Structural sources of robustness in biochemical reaction 
networks. Science, 327(5971):1389-91, 2010. 

[20] M. Thomson and J. Gunawardena. The rational parameterization theorem for multi- 
site post-translational modification systems. J. Theor. Biol., 261:626-636, 2009. 

[21] M. Thomson and J. Gunawardena. Unlimited multistability in multisite phosphory- 
lation systems. Nature, 460:274-277, 2009. 

[22] A. C. Ventura, J. A. Sepulchre, and S. D. Merajver. A hidden feedback in signaling 
cascades is revealed. PLoS Comput. Biol, 4:el000041, 2008. 



Acknowledgements 

This work was supported by the Lundbeck Foundation, the Leverhulme Foundation and the 
Danish Research Council. EF is supported by a postdoctoral grant "Beatriu de Pinos" from 
the Generalitat de Catalunya and the project MTM2009-14163-C02-01 from the Spanish 
government. Part of this work was done while E.F. and C.W. visited Imperial College 
London in 2011. Neil Bristow is thanked for assistance. 



Simplifying Biochemical Models 



9 



Supplementary Information 

Elisenda Feliu, Carsten Wiuf 

This section provides the proofs and explanations of the facts described in the main 
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1 Proofs of facts 
1.1 Preliminaries 

Reaction networks. General standard background material on reaction networks can 
be found in [26,25]. Here we recapitulate the definitions and properties necessary for our 
work. Consider a set S of n species Si, ... , S n . A reaction network (or simply network) 
consists of a set of reactions 1Z whose elements take the form y — > y' with y = Ya=i a i^i 
and y' = Y^i=ifii,jSi f° r some non-negative integer coefficients aj,/3j > 0. The linear 
combinations y, y' are called complexes and the coefficients are called stoichiometric coeffi- 
cients. Complexes y, y' can be seen as elements of the vector space W 1 with entries given by 
the stoichiometric coefficients. An intermediate Y satisfies that the only complex involving 
Y is Y itself and there is at least one reaction of the form y — > Y and one reaction of the 
form Y — > y . Here y and y' can be other intermediates. An intermediate is thus both a 
species and a complex. 

The molar concentration of species Si at time t is denoted by q = Cj(i). To any 
complex y we associate a monomial c y = Yl™—i C T ■ F° r example, if y = (2, 1, 0, 1), then the 
associated monomial is c v = c\c2C^. In the main text, concentrations are denoted by [Si] 
and the monomial associated to y by [y]. 
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We assume that each reaction y — > y' has an associated positive rate constant k y ^ y i, 

t licit IS . ky yyt is in IR + . The set of reactions together with their associated rate constants 

give rise to a polynomial system of ordinary differential equations (ODEs) taken with 
mass-action kinetics: 

6i= ^2 k y^y' cV (y'i ~ 2/i)> i=l,...,n. (7) 

These ODEs describe the dynamics of the concentrations q in time. The steady states of 
the system are the solutions to a system of polynomial equations in c\ , . . . , c n obtained by 
setting the derivatives of the concentrations to zero: 

= ^ k y^y' cV (yi ~ 2/i)> for alH = 1 . . . , n. (8) 

y^y'eTZ 

It is convenient to treat the rate constants as parameters with unspecified values, that 
is as symbols. For that, let 

Con = {k y ^ y t\y -)■ y G TZ] 

be the set of the symbols. Then the system Q is a system of polynomial equations in 
ci, . . . , c n with coefficients in the field R(Con). 

The dynamics of a reaction network might preserve quantities that remain constant 
over time. If this is the case, the dynamics takes place in a proper invariant subspace of 
R n . Let x ■ x' denote the Euclidian scalar product of two vectors x,x' and IR + the vectors 
with non-negative coordinates. 

Definition 9. The stoichiometric subspace of a reaction network with reactions set TZ is 
the following subspace of K n : 

T = (y'-y\y^y'£K). 

By the definition of the mass-action ODEs, the vector c points along the stoichiometric 
subspace T. The stoichiometric class of a concentration vector c is {c + T} n K, . Two 
steady states c, d are called stoichiometrically compatible if c — d G T. This is equivalent 
to u ■ c = oj ■ d for all oj G T 1 - . 

In other words, if oj = (Ai, . . . , X n ) G T- 1 , then Y17=i ^»^» = ^- This implies that the 
linear combination of concentrations Y27=l ^» c * ^ s independent of time and thus determined 
by the initial concentrations of the system. Such a relation is called a conservation law 
and the value it takes in a stoichiometric class is called a total amount. In particular, any 
steady-state solution of the system preserves the total amounts. 

Graphs. Given a directed graph G we call r a spanning tree of G if r is a directed 
subgraph of G with the same node set as G, and the undirected graph obtained by removing 
orientations from edges in r is connected and acyclic. A spanning tree r is said to be rooted 
at v if v is a node in r, and the unique path from any other node w G r to v is directed 
from w to v. G is strongly connected if for any (unordered) pair of nodes v,w G r there 
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is a directed path from v to w. If G is labeled then any spanning tree r will inherit the 
labelling from G in the obvious way. For any labeled graph G we define 

tt(g)= n a - 

Core and extended models. Consider a core model with species Sc = {Si, . . . , S n }, 
set of reactions TZc and let Cc denote the set of core complexes. An extension model (of 
the core model) has the following form: 

(i) The set of species is Se = Sc U y with y a set of intermediates. Let p be the 
cardinality of y. 

(ii) The set of reactions y — > y' obtained from collapsing the reaction paths in the exten- 
sion model y — > Y\ — > ■ ■ ■ — > Yk — > yf with Y G y and y, y' G Cc equals TZc- 

The set of reactions TZe is divided into four non-overlapping subsets: 

- The reactions that are both in the extended and in the core model, TZcnE = IZe^TZc- 

- The reactions from a core complex to an intermediate, IZc^E = {y — > Y <E TZe\ yG 

c,y £ y}. 

- The reactions from an intermediate to a core complex, IZe^c = {Y — > y G 1Ze\ V G 

c,Y€y}. 

- The reactions between two intermediates, IZe^e = {Y — > Y' G 1Ze\ Y,Y' € y}. 

We assume that the set of species of an extended model is ordered as {Si, Yi, ... , Y p }. 

For simplicity, we let Cj denote the concentration of S, for i = 1, . . . , n and U{ the concentra- 
tion of Yj for i = 1, . . . ,p. The ODEs of the extended model consist of n+p equations. Since 
intermediates do not interact with species Si, the ODE equations do not have monomials 
involving both c* and u* . 

1.2 Proof of Fact 1 

Consider a core model with species Sc = {Si, . . . , S n }, set of reactions TZc an d let Cc 
denote the set of core complexes. Let Tc be the stoichiometric space. Consider an extension 
model with set of species is Se = Sc U 3^ with y = {Yi, . . . , Y p } a set of intermediates, and 
set of reactions TZe- Let Te be the stoichiometric space of the extended model: 

r s = (y' -y'\y^y' en E ). 

For every reaction y — > y' G TZc, there exists a reaction path y — > Y x — > ■ ■ ■ — > Yi k — > y' , 
possibly with empty set of intermediates, such that each reaction belongs to TZe- It follows 
that there is an inclusion 

T c T E (10) 
obtained by setting the coordinates n + \, . . . ,n + p to zero. 
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Let the reaction graph of a network be the graph with the complexes as nodes and an 
(undirected) edge between any two complexes forming a reaction. Let the reaction graph 
of the core model have J components. Then the reaction graph of the extension model also 
has J components. Any reaction in the core model can be realized as a series of reactions 
in the extension model, by assumption. Hence the extension model cannot have more 
than J components. We show that it has precisely J components. Consider intermediates 
Yi x ,..., Yi k such that y — Yi x — . . . — Yi k — y' is a series of reactions (here — is either — > or 
<—) and y, y' belong to different connected components of the core reaction graph. If the 
reactions are all in the same direction then either y — > y' or y' — > y is in the core model 
and hence y, y 1 belong to the same connected component of the core reaction graph. If the 
reactions are in different directions, let K . be the first intermediate such that — ¥ Y%. or 
<— Yi . — >. By hypothesis, there exists a reaction path Y{. —>-•■■ — >• y" or y" —>■••• — >• Y^. 
respectively. Then, either y — > y" and y' — >• y" or the reverse reactions are core reactions 
and hence y, y' belong to the same connected component. 

The statement of Fact Q] is: 

Fact 1. The conservation laws in the core model are in one-to-one correspondence with 
the conservation laws in the extension model. The correspondence is obtained by adding 
the same linear combination of the concentrations of the intermediates to the conservation 
laws of the core model. 

Fact [T] will follow from the lemmas below. 

Lemma 11. Assume that the reaction graph of the core model has J connected components 
(which we order) and for j = 1, . . . , J, select a complex y 3 in each component. Let oj = 
(cji, . . . ,oj n ) G and define a,j = oj ■ y K Define a vector oj G M. n+P such that 



\Ui for i = 1, . . . ,n, 

I aj, if Yi is in the j-th component and i = n + 1, . . . , n + p. 
We have 

(i) ul g r£. 

(ii) If oj 1 , . . . ,u d form a basis ofTjj then w 1 , . . . ,uj d form a basis ofT^. 

Proof. First of all, we check that a,- is independent of the choice of y 3 . Fix a component 
Cj of the reaction graph of the core model. For any reaction y — > y' in Cj, we have 
w " (y' ~ y) = and hence u • y' = to • y. Since Cj is connected, a,j is independent of the 
choice of y 3 . Note that if y is a core complex, then to ■ y = to ■ y. 

To show (i), we need to show that to ■ (y' — y) = for all y — > y' G He- Since to G Tq, 
the equality clearly holds if y — >• y' G TZcnE- Consider y — > Yi G 1Zc->e- If it belongs to 
the j-th component, then we have to ■ Yi = dj = to ■ y = u ■ y. Therefore, to ■ (Y^ — y) = 0. 
Similarly we check that to is orthogonal to all reactions in IZe^c an d 7Ze->e- This proves 
(i)- 

To prove (ii) note that if oj 1 , . . . , io d are linearly independent then so are oj 1 , . . . ,oj d . 



Further, by the inclusion (10), dim(rc) < dim(r^). Consequently, 

d = dim(r^) > dim(r^) > d, 
from where it follows that dim(r^) = d and hence to 1 , . . . , oj d is a basis of Tg. □ 
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Lemma 12. For uj = (uji, . . . , u; n+p ) G Pg, define uj n = (uj\, . . . ,w n ). We have 
(nj If uj 1 , . . . ,uj d form a basis of then uj llT , . . . , tJ d7r /orm a basis of T^,. 



Proof. Any reaction y — > y' G T^c satisfies y' — y & Te under the inclusion (10). Hence 
% " (y ~~ 2/') = 0- Since any core complex y has coordinates n + 1, . . . , n + p equal to zero, 
uj ■ y = uF ■ y. This proves statement (i). 

To prove (ii) we use that dim(r^) = dim(r^) (see previous proof). Let u € and 
consider to G as defined in Lemma 11 . Since uj 1 , . . . , uj d form a basis of Lg, we have 



2 = Xiuj 1 + . . . + X d uj d 
for some Aj. Since a; = a; 71 ", be projecting onto the first n coordinates we obtain 

u = Aiw l7r + . . . + X d u> dn . 
Therefore, uj 1w , ...,UJ d7T generate T^, and hence they form a basis. □ 

Note that the constructions of the two lemmas above give the desired correspondence 
between conservation laws since for all uj G T^, we have uj = uj n and for all uj G we have 
uj = (uj) n . 

Remark 13. The results in this subsection show that core and extension models have the 
same deficiency [25]. The deficiency of a network is defined as the number of complexes 
minus the dimension of the stoichiometric space minus the number of connected components 
of the reaction graph. We have proved that the core and any extension model have reaction 
graphs with the same number of connected components, and that both the dimension of 
the stoichiometric space and number of complexes of an extension model increase by the 
number of intermediates. As a consequence, the deficiency remains invariant. 

1.3 Proof of Fact H 

The proof of Fact [2] relies on ideas introduced in [20] and developed generally in [7] . Let 
us recall its statement with the notation introduced above: 

Fact 2. The system of equations Ui = for all intermediates Yi in the system can be 
solved in terms of the core species and is expressed at steady state as a linear sum 
U{ = Hi,yC y ■ A monomial c y appears in the expression if and only if there is a reaction 
path y —>■...—> Yi involving exclusively intermediates. 

Proof. Let us consider the steady-state equations Ui = for i = 1, . . . ,p corresponding to 
the intermediates. These equations take the form 



°= Y k V^Y 1 C y + Y ky^YiUj-i YI k Y>^y+ Y ky i^ Y i 

y^Yi&lc^E Yj—^Yi&'R-E—fE \Yi^y&l E ^ c Y^Y^TIe^e 

(14) 



Ui 



Simplifying Biochemical Models 



14 



(here, i is fixed and summation is over Yj and y). It follows that equations (14) for 
i = 1, . . . ,p form a system of linear equations in the variables u\ , . . . , u p and coefficients in 
M[ConU{ci, . . . , c n }]. That is, equations (14) for i = 1, . . . ,p form the linear system 

Au + z = (15) 

with u = (m, . . . , Up), and A = {a^}, such that for i ^ j we have 



fcy,.-^ if Yj -> e 72.£_kE 
otherwise, 



and for i = j we have 



-ei-di, with ej= fc r^y fc , ^ = ky^y 



We define z = (z\, . . . , z p ) to be the independent term: 



Zi= k y^Y t c y . 

y^Yi&lc^E 

All coefficients but a^j are positive. Further, a^j € M[Con] while zi G M[ConU{ci, . . . , c n }]. 
The column sums of A are not all zero. Indeed, the sum of the entries in column i is 
Ylj=i a j,i = a j,i -ei-di. Note that for i fixed, 

YI a j,i = Y ky ^ Y i = e *- 

Therefore, we have that 

v 

Y a j,i = ~ d i- ( 16 ) 
j=i 

Since by assumption TZe-^c is n °t em Ptyj di ^ for some z and thus the column sums of 
A are not all zero. 

Consider the labeled directed graph Gy with node set y U {*}. We order the nodes 
such that Yi is the i-th node and * the (p + l)-th node. The graph Gy has the following 
labeled directed edges: 

• Yj — ^> Yj if djj 7^ and i ^ j, 

• Yi — h * if di 0, and 

• * ^ ^ if m ^ 0. 

All labels are in M[ConU{ci, . . . , c n }] and are either zero or polynomials in ConUjci, . . . , c n } 
with positive coefficients. By definition of intermediates, the graph Gy is strongly con- 
nected. Indeed, for every intermediate Yi £ y there is a reaction path Y — > Yj 1 —>■■••—> 
Yj t — > y' with y' ^ y and a reaction path y —> Yj x —)■•••—> Yj ; — > Y for some y y. 
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Therefore, there is a directed path in both directions between each intermediate and * in 
Gy, hence also between any two intermediates. 

Let C = {\i,j} be minus the Laplacian matrix of Gy. If i, j < p, then Xij = a^j. The 
entries of the last row of C are \ p +i,i = di for i < p and the entries of the last column are 
Ai, p +i = Zi for i < p. By the Matrix- Tree theorem [28] we conclude that 

in particular, since the (p + l,p + 1) principal minor of C is exactly A, we have 

a:=(-lf det(A) = (-irC {p+ljP+1) = £ tt(t). (17) 

ree(*) 

Since no spanning tree rooted at * can involve a label Z{, a is in fact a polynomial in 
R[Con]. Since Gy is strongly connected, then there exists at least one spanning tree rooted 
at *, and hence (—l) p det(A) is non-zero in M[ConU{ci, . . . , c n }]. It follows that the system 
Au + z = has a unique solution in E(ConU{ci, . . . , c n }). 

For i = 1, ... ,p, we let <7j be the following polynomial in c\, . . . , c n , 

a t = (-iy +1 £ (p+M) = ^ 

ree(y,) 

which is either zero or has positive coefficients in E[ConU{ci, . . . , c n }]. By Cramer's rule, 
we have 

Ui = <fii{Cl, . . . ,c„) = — — — = — , t = l,..., p. 

(-l) p £(p+l,p+l) p 

Since Gy is strongly connected, there exists at least one spanning tree rooted at Yi, and 
o; L / as a polynomial in M[ConU{ci, . . . , c n }]. 

Since a is a polynomial in R[Con], then Ui = Oijo can be seen as a polynomial in 
K[ci, ... ,Cn\ with coefficients in E(Con). Further, each term o; L can be written as: 

v v 

k=i k=i y^Y k en c ^E 

with afc f j € M[Con]. Specifically, a^j is a sum of terms obtained from the spanning trees 
rooted at Yi containing the edge * — > Y^. Each spanning tree gives a term, namely the 
products of its labels, except the label Zk for the edge * — > Y^. If we define 

v , 
fc=i 

(with k y ->Y k = if the reaction y — > Y^ does not exist) then 
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This proves the first part of the statement. 

To prove the second part, we show that the coefficient \ii^ y can be obtained from a 
graphical procedure. For a fixed core complex y, let Gy be the labeled directed graph with 
node set y U {*} and nodes ordered as above. The graph Gy has the following labeled 
directed edges: 



• Yj — ^> Yi if dij 7^ and i / j, 



Yi%*ifdi^ 0, and 



» Yi if k y -+ Yi / 0. 

That is, Gy and Gy have the same edges and differ only in the label of the edges * Yi, 
i = 1, . . . ,p. Then 

= ^ - (19) 

where O y (-) refers to the spanning trees of rooted at the argument. We have that 
^i,y 7^ if and only if there is a spanning tree rooted at Yi in Gy. Equivalently, if and only 
if there exists a reaction path from y (that is, *) to Yi. □ 

1.4 Proof of Fact M 

Let us recall the statement of Fact [3l 

Fact 3. After substituting the expressions Ui = Yl y Mi,2/ c2/ *nfo the ODEs for di of the 
extension model, a mass-action system for the core model is obtained with rate constants 
that are derived from the reaction paths connecting the complexes in the extension model. 

Proof. The system of equations that describes the mass-action kinetics of the core model 
for some constants t y ^ y i is: 

(H= y W c %i-^)- ( 2 °) 

The ODE corresponding to Cj, % = 1, . . . , n, of the extension model taken with mass-action 
kinetics is 

v v 
Ci= Y k y^y'C y (y'i-yi) + Y YI ^Y^y'Ujy'i-^2 Y k y ^ Yj c y yi. 
y-ty'eR-CnE 3=^Yj^y'e_n E ^ c j=l y^Yjellc^E 



Using (18), we obtain 



Ci= Y k y ^ ylC y {y'i-yi) +Y Y k Yj-*y' Y Vja^y'i-Y Y k y ^ Y] c y y l 
y^y'&lcnE j=lYj-+y'eK E ^c y&C c j=l y-+YjeKc^E 

(21) 
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We want to see that this expression can be written in the form of (20) for some choice of 
constants t y ^ y > expressed in terms of A;*. Let 

p v v 

ky^-y' = ^ > ky._ >y' l^j,yi = ^ ] ^ ] ky^ y iC y yi, Bi = ^ ] ^ ] k y ^,y.C y yi. 

i=i 3=ly,y'eC c j=i 2/->K 3 e&c->£ 



where k y ^ y / might be zero if ky._> y i = or = 0. Then (21) can be written as: 

ci = k y^y' cV (y'i -Vi)+ ^2 k y^y' cy (y'i ~ f») + ^ ~ 

y^y'&llcnE y,y'&C c 



Assume that for all fixed i we have = Bi (proven below). Then (21) reduces to 

Ci= ky^ y tc y (y' i -y i )+ ^ k y ^ y > c y (y- - y { ) . (22) 

y^y'&TicnE y,y'&C c 

Let us see that k y ^ y i 7^ if and only if there is a reaction path from y to y' involving 
exclusively intermediates. If (j,j t y 7^ then there is a spanning tree in Gy rooted at Yj. 
In particular, there is a reaction path from y to Yj involving intermediates. If further 
ky^yi 7^ then there is a reaction Yj — > y' which all together give a reaction path y to y' . 
By hypothesis, the reaction y — > y' is in the core model. 

Reciprocally any reaction y — ¥ y' in the core model appears in at least one reaction path 
y — )■ Yjj —)••••—)• 3^, — > y' , potentially without intermediates. If the reaction itself is not 
in the extended model, then ky ik ^y' 7^ and there is a directed path from * to in the 

graph Gy. Since Gy is strongly connected by hypothesis, any such path can be extended 
to a spanning tree of Gy rooted at Y^. It follows that for all reactions y — > y' G IZc \ T^E 
there exists an index k for which ^k, y ky k ^yi 7^ 0. 
Consequently, (22) can be written as 

c% = ^ ] ( k y-^y' k y -> y >)c^ \y^ — yi) 
y^y'ellc 

(with ky^yi = if the reaction y — > y' is not in the extended model). Therefore, by defining 

p 

ty-+y' := k y->y' k y^y' = k y^y' ^ ] ky.^. y i f-tj^y (23) 

3=1 

a mass-action system of the core model is obtained. 

It remains to show that for fixed i we have Ai = B{. It is sufficient to show that for 
fixed y £ Cc with yi 7^ 0, we have 

p p 

3=1 Yj^y'&lE^c 3=1 

where in the right-hand side we allow k y ^y. = if the reaction does not exist. Consider the 
graph Gy defined above. Recall that dj = ^2y ] ^ ye -ji E ^ c ky-^ y and fij }V = ^f-. Therefore, 
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we have to show that for a fixed y S Cq with yi ^ we have 

p p 

Y d 3 a i>y = Y k v^Yi a y ( 24 ) 

j=l 3=1 

Consider the set G\ of all possible subgraphs of Gy which are the union of a spanning 
tree rooted at * and an edge from * to some Yj £ y, and the set G2 of all possible subgraphs 
of Gy which are the union of a spanning tree rooted at some Yj £ y and an edge from Yj 
to *. Observe that we can rewrite ( |24[ ) as 

Y 7F ( T ) = Y 

t£Gi t£G 2 



and so showing that (24) holds reduces to showing that G\ = Gi- 

Let r G G\. There is a single cycle in r, containing at least the nodes * and some node 
Yj, to which the unique outward edge from * points. Along this cycle there is a unique 
inward edge to *, with label d m 7^ for some m. Note that there is a directed path from 
every node in r to *. The directed path from a node w to * either passes through the node 
Y m , or it does not. In the former case, the directed path from w to Y m is preserved if we 
remove the edge from Y m to *. In the latter case, the path from w to * is unaffected if we 
remove the edge from Y m to *, and we can extend this path to (via the edge from * to 
and (if Y^ / Y m ) hence to Y m (via edges which comprise part of the cycle in r). We 
also know that the edge from Y m to * is part of the unique cycle which r contains. Thus 
removing this edge yields a spanning tree of the same node set, but rooted at Y m . Since we 
know that d m 7^ 0, we can add this edge back in to see that r 6 G2. This shows Gi C Gi. 

The proof that G2 Q G\ is analogous, with the roles of * and Yjj reversed. □ 

1.5 Proof of Fact H 

We use the notation introduced in the previous sections. Consider a core model with 
species set Sc and set of reactions TZc ■ Consider an extension model with species set 
<Se = Sc U y with y the set of intermediates, and reaction set TZe- Rate constants t y ^ y i 
of the core model are realizable in the extension model if there exist rate constants k y ^ y > 
in the extension model such that 

p 

ty-+y' = ky^-y' + ^ ] ky^y ' fJ-j,y, (25) 
3=1 

which is the relationship established between parameters in the core and extension model 



in equation ( 23 ) . 

A steady state is said to be non-degenerate if the Jacobian of the ODE system at the 
steady state is non-singular over the stoichiometric space. 
Let us recall the statement of Fact [H 

Fact 4. // the core model has N non- degenerate steady states for some rate constants and 
conserved amounts, then any extension model that realizes the rate constants has at least N 
corresponding non-degenerate steady states for some rate constants and conserved amounts. 
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These can be chosen such that the correspondence preserves unstable steady states with at 
least one eigenvalue with non-zero real part and asymptotical stability for hyperbolic steady 
states. Reciprocally, if the canonical model of a steady-state class can have a maximum of 
N steady states for some rate constants and conserved amounts, then any extension model 
in that class, or in any smaller class has at most N steady states. 

The fact follows from the series of propositions and lemmas below. 

Proposition 26. Consider a core model with species set Sc and set of reactions He- 
Consider an extension model with species set Se = Sc U y with y the set of intermediates, 
and reaction set He- Assume that: 

(i) For some choice of rate constants r = {t y ^. y >}, y — )■ y' G Tic, the core model has 
N > 1 distinct non- degenerate positive steady states in the same stoichiometric class. 

(ii) There exist rate constants k = {k y -> y '} for the extension model that realize t, that is, 
rate constants such that 

v 

i=i 

Then, there exists a choice of rate constants for the extension model that realize r for which 
there are N distinct non- degenerate positive steady states in the same stoichiometric class. 

Proof. We will first rewrite the steady-state equations for the core model and the extension 
model in a way suitable for our purpose. Secondly we show that if the core model has N 
non-degenerate positive steady states in the same stoichiometric class then so does the 
extension model. Let d = dim(r^) = dim(r^) (Fact fl|. We assume that the extension 
model has p intermediates and that the species set Se is ordered as Si, . . . , S n , Y\, . . . , Y p 
where Sc = {Si, . . . , S n } and y = {Yi, . . . , Y p }. We let q denote the concentration of Si 
and Ui the concentration of Y^. 

Consider the core model, a concentration vector c G M" and rate constants r = {t y ^ y '}. 
The steady-state equations are given by 

9t(c) := ^2 ty^y'iv' ~ y)° y = °> 
y^y'eKc 

together with the equations for the conservation laws for a given set of total amounts 
Ti, . . . , T,i. We follow [6] and choose a reduced basis for F c , that is, a basis {ui 1 , . . . , uo d } 
with oj l = (\\, . . . , \ l n ) such that A- = 1 and A* = 0,j / i. Such a basis always exists, 
potentially by reordering the set of species Sc [6] . The system of equations to be solved 
can then be rephrased as 

g T (c) = 0, where g T (c) = (w 1 • c - Ti, . . . , uj d ■ c - T d , g T)d+ i(c), . . .,g T ,n{c)) 

(see [6]). Thus, two vectors c, d G Wt are steady states of the core model, for the rate 
constants r, in the same stoichiometric class if and only if g T (c) = g T {c') = for some 
choice of Ti, . . . , T&. 
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Similarly, consider the extension model, a concentration vector (c, u) € R™ +p , and rate 
constants k = {k y -> y i}. The steady-state equations are given by 

o = f K ( c , u)= k v^y - y) c f ■■■■■ c n< +1 ■■■■■ u p n+p > 

together with the equations for the conservation laws for a given set of total amounts 
T%,..., Td- The conservation laws are related to the conservation laws of the core model 
by Lemma 11 and we use the notation introduced there. It follows that if {uj 1 , . . . ,ui d } is 
a reduced basis for then {to 1 , . . . , ui d } is a reduced basis for F^, and that the system of 
equations to be solved can be stated as 

f K (c, u) = 0, where 

f K {c,u) = {uj 1 ■ (c,u) -T 1 ,...,uj d ■ (c,u) - T d , f Kjd+1 (c,u), . . . , f K , n+p (c,u)). (27) 

Since d < n, the last p components of f K (c, u) are the steady-state equations corresponding 
to ii = 0. Note that 



p 

uj 1 • (c, it) - Ti = oj % ■ c + ^ w l n+j Uj - Ti, 

i=i 

i = l,...,d, where w l n+ j is the (n + j)-th coordinate of co l as defined in Lemma |TTJ 
Two vectors (c,u),(c',u') £ are steady states of the extension model in the same 

stoichiometric class for the rate constants k if and only if f K (c, u) = f K (c', u') = for some 
choice of Ti, . . . , T d . 

We will reformulate the equation f K (c,u) = to obtain a system of equations that is 
closely related to the equation g T (c) = 0. First recall that at steady state Ui = Y^ y ec c A^s/ ' 



y 



(Fact[2j). In equation (27) we will replace the functions / Kj j(c, u), i > n, by the functions 
f K! i(c,u) = Ui — Y^y^Cc ^i,y cy ' ■> i > n 7 an d further replace the variables Uj, j > n, by 

Hyacc ^hv cV in fK,i( c ' u )i for a11 i - n - 

Formally, we proceed in the following way. Let I r denote the identity matrix of order 

r. Note that the function f K (c,u) is linear in u and can be written in block form as 

/«(*«)=( A ) tt+ (x 

where M is an n x p matrix with entries in R[Con], v a vector of length n with components 
in IR[Con, c, Ti, . . . , Td] and A, z are given in the proof of Fact[2j that is, from equation ( 15 ), 
we have that 

(f K ,n+i{c,u), ■ . ■ , f K , n +p{c, u)) = Au + z. 

The p x p matrix A has entries in M[Con] and is invertible in M(Con). The vector z has 
length p and depends on c and Con. Let A" 1 be the inverse of A in IR(Con) . By Fact 
pi the solution to Au + z = is given by itj = — (A~ 1 z)i = Y^ y Hi,yC v ■ Let B be the 
{n+ p) x (n + p) matrix defined in block form by 

p I n -MA' 1 

L { ,1 1 
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This matrix is invertible in R(Con). Then, the function f K (c,u) defined by 

f K (c,u) := Bf K (c,u) 



(28) 



fulfills 
Indeed, 



f K l ci, . . . , c n , J2y&c c Mi>2/ c2/ ' • • • ' Sj/eCc ^p<v c% 



Uj 



B f K (c, u) 



-MA' 1 
A- 1 



Mu + v 
Au + z 



i = 1, . . . , n, 

i = n + 1, . . . , n + p. 



v - MA~ x z 
u + A~ l z 



and the claim follows from the equality — (^4 1 z)i = ^2* y Hi,yC y . 

Note that / K ,i(c, u), i < n, does not depend on u. Further, solving f K (c,u) = is 
equivalent to solving f K (c,u) = 0. Equation (28) ensures that the determinant of the 
Jacobian of f K evaluated at (c, u) is non-zero if and only if the determinant of the Jacobian 
of f K evaluated at (c, u) is non-zero. Consequently to study non-degenerate steady states of 
the extension model we can study zeros of f K (c,u) for which the Jacobian is non-singular. 
This is what we do next. 

Assume that the core model has N positive non-degenerate steady states, d G W\, 
i = 1, . . . , N, in the same stoichiometric class for some rate constants r = {t y ^ y /}, y — > 
y' G "R-c- Let Ti,...,T^ be the total amounts defining the stoichiometric class for the 
reduced basis {a; 1 , . . . , ui d }. 

Let k = {ky^yi} be rate constants for the extension model (23) such that 



Yj->y'l J <j,y 



for all reactions y — )■ y' in the core model TZc (which exist by assumption) 
construction and using Fact [3] we have 



Then by 



/«,i(c,«) 



9r,i (c) + Ej=i <+j E y Vj,yC y i 
5r,i(c), i 



l,...,d 

d + 1, . . . , n. 



Let 6 € be a positive constant. Define a new set of rate constants kP = {ky_^ y ,} by 



ki 



ky^y'/9 if y ->■ y 1 G 11e->e or 11e->c and k e v _^ v , = k y ^ y > otherwise. Let t e y ^ yl and 



f i j v correspond to t y ^ y i and fij^ y , respectively, obtained with the rate constants k using 
(|23j) and §L9§. Then 



^jtV — and ^y^y 1 — l y—>y'' 

The function f®(c,u) for the rate constants kP takes the form 



tit 



igrA^ + oiEU^n+j^yW") 



9r,i( C ) 

Ui-0{EyeC C 



^i,yC y ) 



i = l,...,d, 
i = d + 1, . . . , n, 
i = n + 1, . . . , n + p. 
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We observe that the Jacobian of at (c, it), J[c,u)Uk)i takes the block form 

r / Jc(ffr) + A(*) 



where "(*)" indicates some matrix that we are not concerned with knowing the exact form 
of. 

By continuity, the function f^(c, u) is well defined for all 9 E M. That is, there is a well 
defined and differentiable function 

K" x I p x 1 ^> 

(c,u,0) ^ F«(c,«,e):=g(c,«). 

For = 0, the vectors (c l , 0), 1 < i < JV, are non- negative steady states in the stoichiometric 
class of the extension model defined by the total amounts T\, . . . , That is F K (c l , 0, 0) = 
for all i. The Jacobian of f®(c,u) has the matrix in block form 

J(cu)(U-[ o Ip 

Since the Jacobian matrices of g T evaluated at c l , 1 < i < N, are by assumption non- 
singular, the Jacobian matrices of f®(c,u) evaluated at (c\0) are non-singular. Therefore, 
the Implicit Function Theorem applied to F K at the point (c l ,0, 0) guarantees that there 
exists an interval Ii = (—fa, fa), fa > 0, and an open neighborhood U{ of (c*,0) such that 
for all 9 £ Ii there is a steady state (c l (6),u l {0)) £ Ui in the stoichiometric class defined by 
T%, . . . ,T<i and with (c*(0), u l (0)) = (c*, 0). By making fa sufficiently small, the interval Ii 
can be chosen such that c l {9) is positive (i.e. £/» C WL x M p ) and the Jacobian of f®(c, u) 
evaluated at (c 1 (9) , u l (9)) is non-singular for all 9 £ Ii. Restrict Ii to the positive part, 
if = [0,fa). Since c l {9) is positive if follows from the definition of f®(c,u) that u l {9) is 
positive for all 9 E if . Hence (c* (9) , u l (9)) is a positive non-degenerate steady state in the 
stoichiometric class defined by the total amounts T\ , . . . , for all 9 G if . Since c % ^ c J , 
for all i / j, then by choosing 4> l small enough we are guaranteed that r\f =l Ui = 0. 

With these data, let (f) = mm(^fa\l < i < N). Then, for all 9 E (0, fa) the rate constants 
k s = {ky_iyi}, y —5- y' E TZe, fulfill that the extended model has positive distinct non- 
degenerate steady states (c 1 (9) , u l (9)) , 1 < i < N, in the stoichiometric class defined by 
T\, . . . , Td- This concludes the proof. □ 

Remark 29. A steady state in the core model has always a corresponding steady state in 
the extension model for any choice of matching rate constants, k and r. It follows from the 
following: a steady state c in the core model always defines a steady state concentration u 
for the intermediates. By construction (c, u) is a steady state. If there are steady states 
in the core model in some stoichiometric class for some rate constants then we are however 
not guaranteed that N corresponding steady states in the extension model are in the same 
stoichiometric class. If the stoichiometric space of the core model has full dimension then 
the stoichiometric space of the extension model has full dimension (Fact[l]). Consequently 
the Af steady states are always in the same stoichiometric class. 
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Lemma 30. Let A be the px p matrix in equation (15). Then all p eigenvalues of A have 



negative real part, that is, if X is an eigenvalue of A then Re(A) < 0. 

Proof. We will need the following fact (*): A Metzler matrix M is a square matrix with all 
off-diagonal entries non-negative. If M is a Metzler matrix then exp(M) is a matrix with 
non-negative entries. If M is a Laplacian matrix, then — M is a Metzler matrix, exp(M) 
is a matrix with non-negative entries and all column sums equal to one. This result and 
the Perron- Frobenius theorem used later in the proof can be found in [23]. The argument 
we give holds generally for Metzler matrices with non-positive column sums, but we have 
not been able to find a reference to it in the literature. 

Equation (17) shows that (— l) p det(^4) is a non-zero polynomial in M[Con] with positive 



coefficients. Hence zero cannot be an eigenvalue of A for any choice of rate constants. By 
definition, A is a Metzler matrix, and thus B = exp(A) has non-negative entries. We 
extend A to a (p + 1) x (p + 1) matrix 



A 



A 



(di)i 







where 0„ is the p-dimensional column vector with entries and di are defined in the proof 
of Fact|2l By (16), —A is a Laplacian, and hence B = exjp(A) has non-negative entries and 



all column sums are equal to one. The matrix B takes the form, 



B 



B P 
D 1 



where D is a 1 x p matrix. It follows that the column sums of B are less than or equal to 
one. 

An eigenvalue \x of B is related to an eigenvalue A = Ai + i\ 2 of A by \i = exp(A). 
Assume first that A is irreducible. Hence also B is irreducible. It follows from the Perron- 
Frobenius theorem that all eigenvalues [i of B fulfill < r < 1 (the maximal column sum) 
for some real number r and that [i = r is an eigenvalue. Since A = is not an eigenvalue of 
A, then necessarily r < 1. Hence for all eigenvalues \x = exp(A) of B, we have e Al = |/x| < 1 
and hence Ai = Re (A) < for all eigenvalues A of A. 

If A is not irreducible then A can be written in the following form, potentially after 
reordering the intermediates, 



A 



( M 
o 



V o 



A 2 





\ 



A k J 



where k > 1 and Ai, . . . , Ak are irreducible square matrices. Each Aj fulfills the same 
properties as A above, that is, Aj is a Metzler matrix with non-positive column sums and 
with at least one negative column sum. The latter follows from the following. Let 3^' 
denote the set of intermediates corresponding to the rows of Aj . and let j\ , . . . , jt be the 



corresponding ordered row indices. Using (16) and the definitions above it, and the block 
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diagonal form of A, we have that the column sums of Aj are given by 

it p 

By definition of intermediate, there exists at least an intermediate lj in yj and a reaction 
li — > X with X an intermediate not in or a core complex. As a consequence, there exists 
an index % such that dj 7^ or aj j 7^ for some I < j\. Therefore the column sums of Aj 
are not all zero. Since the eigenvalues of A agree with the eigenvalues of Aj, j = 1, . . . , k, 
the lemma follows from considering each irreducible matrix Aj by itself. □ 

Remark 31. It follows from [6, Remark 7.8] that for a non-degenerate state state, the 
eigenvalues of the corresponding Jacobian matrix can be ordered such that Aj = for 
i = 1, . . . , d and Aj 7^ for i = d+1, . . . , n, where n—d is the dimension of the stoichiometric 
space. 



Proposition 32. Assume as in Proposition 26 Let c = g T (c) be the ODEs describing 
the core model and (c,u) = f K (c,u) the ODEs describing the extension model, where k 
realizes r. Let Aj, j = l,...,n, be the eigenvalues of the Jacobian of g T evaluated at a 
non- degenerate positive steady state d , ordered such that Aj = for i = 1, . . . , d and Aj 7^ 
for i = d + l,...,n. Further, let a 1 , i = 1, . . . ,p, be the eigenvalues of the matrix A in 
equation (15). 

Then n can be chosen such that the extension model has N non- degenerate positive 
steady states in the same stoichiometric class, each corresponding to one of the N steady 
states of the core model. Further, the following holds. Let Vj, j = l,...,n + p, be the 
eigenvalues of the Jacobian of f K evaluated at the steady state (c*,u*) corresponding to d . 
Appropriately ordered the eigenvalues fulfill: 

(i) Vi = 0, i = 1, . . . ,d. 

(ii) If sign(Re(Aj)) 7^ then sign(Re(^j)) = sign(Re(Aj)), i = d + 1, . . . , n. 
(Hi) sign(Re(i/j)) = sign(Re(aj_ n )) < 0, i = n + 1, . . . , n + p. 
Consequently: 

(iv) If a steady state in the core model is unstable and Re(Aj) > 0, for some i = d + 
1, . . . , n, then the corresponding steady state in the extension model is unstable. 

(v) If a steady state in the core model is hyperbolic, that is, Re(Aj) 7^ for i = d+1, . . . , n, 
then the corresponding steady state in the extension model is hyperbolic 

(vi) If a hyperbolic steady state in the core model is asymptotically stable then the corre- 
sponding steady state in the extension model is hyperbolic and asymptotically stable. 



Proof. We use the notation introduced in the proof of Proposition 26 and proceed as in 
the proof of that proposition. The function / K (c, u) is linear in u and can be written in 
block form as 

f K (c,u) = ( A ju + 
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where M' is an n xp matrix with entries in M[Con], v' a vector of length n with components 



in M[Con, c, T\, . . . , Tj\ and A, z are given in the proof of Fact [2j equation ( 15 ). The px p 
matrix A has entries in M[Con] and is invertible in M(Con) with inverse A^^^. The vector 
z has length p and depends on c and Con. Let B' be the (n + p) x (n + p) matrix defined 
in block form by 

B = { A^ 
This matrix is invertible in M(Con) with inverse 

_ ( I n M' 



A 

Then, the function f K (c,u) defined by 

f K {c,u) := B'f K (c,u) (33) 

fulfills 



fn,i{c,u) 



Ik ^Cl, . . . , C n , XyyeCc Ml,2/ C ^> ■ ' ' ! X/j/SCc foiV^J 2 — 1, . . . , n, 

Ui-Y.ytCcV" 1 ^ i = n + l,...,n + p. 



Let k = {/cy^y'} be rate constants that realize t (which exist by assumption). Then 
by construction and using Fact [3] we have 

f K ,i(c,u) =g T ,i(c), i=l,...,n. 

Let £ M+ be a positive constant. Define the rate constants k? = {ky_^ y ,} identically 



to how we did in the proof of Proposition 26 Then the function f^(c,u) for the rate 
constants kP takes the form 



fL(c,u) 



9r,i(c) i = l,...,n, 

Ui - 0{Yjy^c o ^y cV ) i = n + l,-.-,n + p, 



and we observe that the Jacobian of j 6 K at (c,u), </( c ,u)(/k)> is a block matrix 

T ( Jc ^t) 

where "(*)" indicates some matrix that we are not concerned with knowing the exact form 
of. It follows that the Jacobian of f% at (c, u), J( c ,u)(fk) ^ s 

T (f e, o/-i ( U9r) \ _ / J c (g T )-M'9(*) M> 
J(c,u){J K ) o y ^ J y -A6(*) A 

and that the characteristic polynomial X(c,u) 

(x) fulfills 

X{c,u){x) = det(J( C)U )(/f) - xl n+p ) = det(J c (g T ) - xl n ) det(^4 - xl p ) + 9h K (x), 



Simplifying Biochemical Models 



26 



where h K is a polynomial whose precise form we are not interested in. 

We now assume that the core model has N > 1 non-degenerate positive steady states 



in the same stoichiometric class for r. Proposition 26 guarantees that there exists cf> G M+, 
such that for 9 £ (0,4>) the extension model with k° has N non-degenerate positive steady 
states in the same stoichiometric class. Let the steady states in the core model be c l , 
i = 1, . . . , N, with corresponding steady states in the extension model being (c 1 (9) , u l (0)) , 
i = 1, . . . , N (which vary continuously in 9). 



The dimension of the stoichiometric space of the extension model is n+p—d (Lemma 11 ) 



hence there are d zero eigenvalues (Remark 31 ). This proves (i) by ordering the eigenvalues 



appropriately. The characteristic polynomial x is continuous in (x, 9) (it is a polynomial in 
both variables for a fixed steady state and (c(9),u(9)) is a continuous function). Hence the 
zeros of x vary continuously in 9, for 9 small enough. By choosing <p potentially smaller, by 
continuity we are guaranteed that if the real part of an eigenvalue of J c (<?t) is non-zero then 
so is the real part of the corresponding eigenvalue of J c {e),u{6){ft)- Since there is a finite 
number of eigenvalues for the TV" steady states, <j) can be chosen such that the statement is 
true for all eigenvalues. Hence (ii) is true. Similarly </> can be chosen so small that (iii) is 



true. The negativity of Re(fj) for % = n + 1, . . . , n + p, follows from Lemma 30 

Assume now 9 is chosen such that (i)-(iii) are true. Consider an unstable steady state 
for the core model with Re(Aj) > for some i. Then also Re(z^) > according to (ii). 
Positivity of the real part of an eigenvalue implies that the steady state us unstable [27], 
hence the steady state in the extension model is unstable. It proves (iv). A steady state 
is hyperbolic if all eigenvalues have non-zero real part [27]. Then (v) follows from (ii) and 
(iii). A hyperbolic steady state is asymptotically stable if and only if all eigenvalues have 
negative real parts [27]. It follows that if a steady state in the core model is asymptotically 
stable then Re(A«) < for all i = d + 1, . . . , n. According to (ii) we also have Re(i^) < 
for all i = d + 1, . . . , n. Together with (iii) the steady state in the extension model is 
asymptotically stable. It proves (vi). □ 

2 Realization of rate constants 
2.1 Canonical models 

Consider a core model with species set Sc and set of reactions IZc- Consider a canonical 
extension model with dead-end at some core complex y* . That is, the extension model has 
set of species Se = 5cU{Y} (Y is an intermediate), and set of reactions TZe = 7£cU{y* — > 
Y,Y^ y*}. 

Consider some choice of rate constants r = {t y ^ y i}, y — > y' G TZc in the core model 
and total amounts T\, . . . , corresponding to some choice of basis of Yq. We prove here 
that there exist rate constants k = {ky^yi} for the extended model realizing r, that is, 



such that equation ( 25 ) holds 



h-^y 1 ~ ky^y' + 2_. ^Yj^y'^j,y 
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In this case, there is only one intermediate and we have 

ky*— -fY 



fJ>Y,y* 



k 



Hence, for all reactions y — >■ y' in TZc, we have 

ty— >y' — ky^ryi 

and realization parameters obviously exist. 

Further, any conservation law in the extended model that is not a conservation law in 
the core model takes the form 

uj = u) + aY 

for some constant a. Written as an equation in the core species, we have 

~ , ky*^.Y y* 

UJ = UJ + drjf C y . 

ky-^y* 

We note that by varying the two rate constants k y *-j,Yi ky-ty* > the coefficient of c y 
takes any desired non-zero value (if o / 0). 

2.2 Non-realizable constants 

Consider the core model with reactions: 




Consider the following extension model: 




Then we claim that this extension model cannot realize all choices of rate constants of 
the core model. If all rate constants of the core model, t\, . . . ,t$ were realizable, then we 



could find rate constants k\, . . . ,k& such that equation (25) holds, that is 



h 



&3 + &4 + h ' 

fc2&3 
&3 + &4 + &5 ' 



*2 



ki k. 



IK4 



k 3 + k^ + fc 5 
k 3 + k^ + k 5 



*3 



fei k. 



1«5 



h + ki + k 5 ' 

k2h 
&3 + k 4 + k 5 ' 



(34) 
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Choose for instance 

ti = 3, t 2 = 4, * 3 = 5, t A = 6, t 5 = 8, t 6 = 15. (35) 



Using ti,tz and (34) we see that k 2 = 1k\. Using t% and t§ we see that k 2 = 3A;i and hence 
system (34) has no positive solution. 

This conclusion can also be derived by noting that the core model has six independent 
parameters, t\, . . . , kf, whereas the extension model has only five, k\,...,ks. 

2.3 Deciding on realizability of rate constants 

In some cases, manual inspection suffices to decide whether an extension model can realize 
all choices of rate constants for the core model. However, it would be desirable to have an 
automated procedure to decide this. 

We give here a necessary criterion that makes use of computational algebra tools, 
namely, Grobner bases. The realizability problem can be stated as follows. Let mc, we be 
the number of reactions in the core and extension models respectively. Consider the map 

J 

{ky^yl} I-)- < ky-^yl + ^ ky^y'Hj,'. 

(assuming that the reaction sets are ordered). Asking for all choices of rate constants in 
the core model to be realizable in the extension model is equivalent to requiring that T is 
a surjective map over the positive orthant M T ? E . A minimal criterion is that tue > mc 



(see Section 2.2 for an example). 

Since /x* are rational functions in the rate constants fc*, T extends to a rational map 
over M. mE and the Zariski closure of the image of T is a real algebraic variety that is defined 
by some ideal / of . . . ,t mc ] [24]. That is, Im(T) = V{I). If / is not the zero ideal, 
then T is not surjective over the positive orthant. Thus, for T to be surjective, a necessary 
condition is that I is the zero ideal. 

The ideal I can be obtained using the Implicitization procedure as described in [24] 
§3. Let fei, . . . , k mE , t\, . . . , t mc be variables corresponding to the rate constants in the 
core and extension models, respectively. Let T = (Ti,...,T mc ) and write the compo- 
nents Ti as a quotient of polynomials in k%,. .. ,k mE : Ti = fijg%. Let J be the ideal of 

k\ , . . . , k mE , t\ , . . . , trac 

] given by 



J — (dlh — fl, ■ ■ ■ ,9m c tm c — fm c , 1 ~ 9l • ■ ■ ■ ' 9rn c 



The last polynomial can be dropped if all gi = I and if some g, t are repeated, we consider 
them only once. Then / is the elimination ideal 



i = JnR[ti, ...,t 



A set of generators of I is given by the polynomials involving t\, . . . , t mc only in the Grobner 
basis of J with the lexicographical order on the order of variables z > k\ > ■ ■ ■ > k mE > 
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t\ > ■ ■ ■ > t mc . Therefore, if such a Grobner basis has no polynomial in t\, . . . , t mc , then 
any choice of rate constants of the core model is realizable in the extension model. 



Example. Consider the example given in subsection 2.2 The map T is 



» 5 



/ci k 



kik^ 



kih 



k 2 k 3 



k 2 k 4 



kok. 



2«-5 



(k±, ks) | y 

\k 3 + k 4 + k 5 ' k 3 + k 4 + k 5 ' k 3 + k 4 + k 5 J k 3 + k 4 + k 5 7 k 3 + k 4 + k 5 : k 3 + k 4 + k 5 

It is clear that this map cannot be surjective over the positive orthant, but in general this 
will not be the case. The ideal J is in this case: 

J = (ti - k\k 3 ,t 2 - k\k^t 3 - kik^,U - k 2 k 3 ,t5 - k 2 kd_,t 6 - k 2 h, 1 - (k 3 + k 4 + k 5 )z). 

Using Maple, we compute the Grobner basis Gj of J with the lexicographic order that 
orders A;*, z larger than i* and obtain: 

Gj ={-t 5 t 3 + t e t 2 , -Ut 3 + t 6 ti, -t 4 t 2 + t 5 ti, -t 5 k 5 + t G k A , -k 5 t 2 + t 3 k A , -t 4 k 5 + t 6 k 3 , 

- Uk^ + t 5 k 3 , -k 5 h + t 3 k 3 , -k 4 ti + k 3 t 2 , -t& + k 2 k 5 , -t 5 + k 2 k 4 , -t 4 + k 2 k 3 , 

- t 3 k 2 + teki, -t 2 k 2 + t$ki, -t\k 2 + tiki, -t 3 + fei&s, — *2 + k\k 4 , -t\ + /ci/c 3 , 

Ztg + Zts + Zt4 — k 2 , Zt% + Zts + ztd_ — kl, — 1 + + + 

It follows that 

I = (~t 5 t 3 + t 6 t 2 , -t 4 t 3 + t 6 ti, -Ut 2 + t 5 tx) ± 

and hence T is not surjective and there exist non-realizable rate constants. Observe that 
the polynomials in / do not vanish when evaluated in the rate constants in ( 35 ) . 

Consider again the core model in subsection 2.2 but now with the extension model 
given by 



2/2 ' fc 2 

The map T is given by 

™S T. m a 





(ki,...,k 8 ) 



fci k 



1*3 



k\ki 



+ 



k^h 



7 



hk 



5«-8 



k 2 k 3 



kok 



2«-4 



+ 



k 6 k 7 



kek 



6^8 



k 3 + k 4 k 3 + &4 k 7 + k & kj + A; 8 fc 3 + /c 4 fc 3 + k A k 7 + k 8 k 7 + fc ; 
The ideal J is in this case: 

J = {ti-kik 3 ,t2-k 1 ki-k 5 k 7 , t 3 -k 5 k 8 , t 4: -k 2 k 3 , t 5 -k 2 kd_-k 6 k 7l t G -k 6 k 8 , (k 3 +k 4 )(k 7 +k 8 )z) . 

We proceed as above and compute the Grobner basis Gj in Maple. In this case, we obtain 
that 

7 = 0, 

indicating that all rate constants might be realizable in this extension model. In fact, in 
this case we find that 

t 3 for t\ k\ . . kd . . k 7 

r = ir> r = ir, t 2 + ts = {t l + u)-^ + {t 3 + t (i ) 1 L , 

t 6 k 6 td_ k 2 k 3 k 8 



which has a positive solution (k\ 



, ks) E M§_ for any positive choice of (ti, . . . , te) E 
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3 Information on the figures 

3.1 Figure 2. Computation of the steady states 

We consider a choice of complexes yi,y2, U3 in Figure 2 that represents a two-site phospho- 
rylation event and compute the maximal number of steady states that the core model and 
the canonical models in Figure 2 can have. Through this section, we are only interested 
in positive steady states and hence, when saying steady state we implicitly mean positive 
steady state. 

Specifically, we consider a substrate S that has two phosphorylation sites, with phospho- 
rylation and dephosphorylation being sequential. We let 5o denote the unphosphorylated 
substrate, Si denote the substrate with the first site phosphorylated and S 2 denote the 
fully phosphorylated form. Phosphorylation reactions are 

S + E^hS! + E S 1 + E%S 2 + E S + E%S 2 + E (36) 

where E is a kinase. By setting y% = So + E, y 2 = S\ + E and y 3 = S 2 + E this model is an 
instance of the core model in the main text, Figure 1A. This model however accumulates 
at steady state all the substrate concentration in S 2 . In order to have a more interesting 
analysis, we add simple dephosphorylation reactions 

S 2 Si % S . (37) 

We want to determine how many steady states can the different canonical representa- 
tives of each class have. If we were simply interested in determining if the system can have 
multiple steady states or not, then we could use one of the several available automatized 
methods (e.g. [5,6]). 



Core model. The ODEs of the core model (Figure 1A) with reactions in (36) and 



(37) are the following: 



[So] = ~h [So] [E] - k 3 [So] [E] + k 5 [5!] , (38) 

[Si] = -k 2 [Si][E] - k 5 [Si] + h[S ][E] + k A [S 2 ], (39) 

[S 2 ] = -k A [S 2 ] + k 2 [Si][E] + h[S ][E], (40) 

[E] = 0. (41) 

This system has two conservation laws: 

Stot = [So] + [Si] + [S 2 ], E tot = [E], 

that is, the concentration of kinase is clearly constant. 

The steady-state equations are obtained by setting the left-hand side of the ODEs to 



zero. Using the steady-state equation derived from (38) we obtain that 



rol (fci + h)E tot 
[oi\ = ^ [ool 
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and then using this expression and MOl) we have 



f k 2 (h + k 3 )Ef ot k 3 E tot \ 

[S2] = v — hh — + ) [o] - 

Using the total amount Stot we obtain that at steady state 

rcl G f-t , hh + k 3 k 4 + k 3 k 5 k 2 (h + k 3 ) 2 

[So] ~ Stot V + i^h Etot + k,h Etot 

Given any positive rate constants k* and positive total amounts E to t and Stot, [So] is 
positive and uniquely determined at steady state by this expression. Further, the steady- 
state value of [Sq] determines the steady-state values of [Si] and [S 2 ] using the expressions 
above. We conclude that the core model has one positive steady state for each choice of 
rate constants and total amounts. 

Extension model 1. We consider the canonical representative of the class in the 
first column of Figure 2. The reactions of the model are those in (36) and (37) together 
with 

S + E %Y, Y % S + E. 
The ODEs of the model are the following: 

[S ] = -h[S ][E] - k 3 [S ][E] - k 6 [S ][E] + fc 5 [5 a ] + k 7 [Y], (42) 

[Si] = -k 2 [S!][E] - k 5 [St] + kt[So][E] + h[S 2 ], (43) 

[S 2 ] = -h[S 2 ] + k 2 [S!][E] + k 3 [S Q ][E], (44) 

[E] = -k e [S ][E] + k 7 [Y], (45) 

[Y] = -k 7 [Y] + k 6 [S ][E]. (46) 

This system has two conservation laws: 

Stot = [So] + [Si] + [S 2 ] + [Y], E t ot = [E] + [Y]. 

Observe that the conservation laws of this system and the conservation laws of the core 
model are in correspondence, as indicated by Fact[T] Isolating [E] from the kinase conser- 
vation law we have 

[E] = E t ot - [Y], 



which is positive provided < [Y] < E to t- From the steady-state equation derived from 

\Sa 



(46) we have that 

k 7 [Y] 



k 6 (E tot -[Y]Y 



Using the steady-state equation corresponding to (42 )+(|46[) and (44) we iteratively obtain 

(fci + fc 3 )fc 7 _ k 3 k 7 Mki + k^ 

[Sl] ~ k 5 k 6 [n [S2] ~ W 1 + k,k 5 k 6 [Y]{Etot ~ [Y]) - 
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Given < [Y] < Etot, all the steady-state expressions above are positive. The value at 
steady state of [Y] is found by imposing the substrate conservation law (Stat) to be fulfilled: 

K6(-Etot - J J k 5 /c 6 fc 4 A; 6 k 4 k 5 k 6 

(47) 

Let us focus on the right-hand side of this expression such that Stat = ^{[Y])- The function 
ip is continuous for < [Y] < E to t and tends to infinity as [Y] approaches E to t- When 
[Y] = we further have tp(0) = 0. It follows that for any given Stot there exists [Y] G 
[0,Etot) such that Stot = <p{[Y]) and hence a positive steady state exists. 

If the function (p is always increasing, then there is exactly one. If it can decrease in 



some part, then there can be more than one. Note that equation (47) can be rewritten 
as a polynomial of degree 3 in [Y] such that the roots in < [Y] < E to t are the positive 
steady states. It follows that there can be at most three positive steady states. 



Each summand in (47) is an increasing function of [Y], except for the last summand. 
Hence, it is not clear whether <p can be decreasing in some interval. The derivative of (p 
with respect to [Y] is 

tf\v}\ 1 i kl ( Etot , fc i + fc 3 , fe 3 k 2 (h + k 3 ) 

<p y = 1 + — -= p-r^ + + — + 7-r {Etot ~ 2 7 

k 6 \(Etot ~ [Y\y fc 5 k 4 k 4 k 5 

This derivative is negative if and only if 

1 , k 7 ( Etot , h + h | k 3 \ k 7 k 2 (h + fc 3 ) 

l + k 6 {(Eto t -[Y]r + h + k A ) < h k A k 5 [2[Yl ht ° t} - 

The left-hand side of the inequality is an increasing function in [Y] that tends to infinity 
as [Y] approaches Etot- The right-hand side of the inequality is a line with positive slope 
that takes a negative value at [Y] = and crosses the x-axis at [Y] = Et t/2- If the line 
intersects the left-hand side curve, then there will be multiple steady states. The left-hand 
side does not depend on k 2 while the slope of the line increases with increasing k 2 . By 
fixing all constants except k 2 and letting k 2 vary arbitrarily, the two curves must meet. 
Except if they meet tangently, the two curves will cross in two points, between which (p 
decreases. In this case, <p increases initially, decreases for some interval, and increases to 
infinity afterwards. It follows that there are values of Stot for which the system has three 
steady states. 

Specific rate constants for which the system has three positive steady states are: 

ki = l,i^2, k 2 = 2, E tot = 10, Stot = 100. (48) 

The three steady states correspond to [Y] = 6.5 — \/IT, 8, 6.5 + \/TT. 

Extension model 2. We consider the canonical representative of the extension 



model class in the second column of Figure 2. The reactions of the model are those in (36) 



and (37) together with 



So + E%Y u Y 1 %S + E, S X +E%Y 2 , Y 2 %S 1 + E. 
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This model can be seen as an extension model of the canonical representative in column 1 
(which is taken as the core model). Canonical models always realize parameters of the core 
model (see Section [2]). Therefore, since the canonical representative in column 1 admits 
multiple steady states, then so does the canonical representative in the second column of 
Figure 2 and it has at least 3 steady states. 

To show that it has at most 3 steady states we proceed as above. We consider the ODE 
system and we iteratively eliminate variables to obtain that at steady state 

k 7 k 9 E tot k 5 k 7 [Si](k 9 + k 8 [Si]) 

h = T. L , T. L ro 1 , L Z rc t l b °\ ~ 



k 7 k 9 + k 6 k 9 [S ] + k 7 k 8 [Sx] ' kg((ki + k 3 )k 7 E tot - k 5 k e [Si\) ' 

and 

*- = W + PJ + ((^ + ^)[ai + (^ + ^)w)w. 

By writing the expression above as a polynomial in [Si], we obtain a polynomial of degree 
3 and hence at most three positive steady states can occur. 

Extension model 3. We consider the canonical representative of the class of exten- 
sion models in the third column of Figure 2. The reactions of the model are those in (36) 
and (37) together with 

Si + E%Y U Yi%Si + E, S 2 + E%Y 2 , Y 2 %S 2 + E. 

The ODEs of the model are the following: 

[So] = -h [S ] [E] - k 3 [So] [E] + k 5 [Si] , (49) 

[Si] = -k 2 [Si][E] - k 5 [Si] - h[Si][E] + ki[S ][E] + k 4 [S 2 ] + k 7 [Yi], (50) 

[S 2 ] = -k A [S 2 ] + k 2 [Si][E] - k 8 [S 2 ][E] + k 3 [S ][E] + k 9 [Y 2 ], (51) 

[E] = -k 6 [Si][E] + k 7 [Yi] - k 8 [S 2 ][E] + k 9 [Y 2 ], (52) 

[Yi] = k 6 [S 1 ][E]-k 7 [Yi], (53) 

[Y 2 ] = k 8 [S 2 ][E] - k 9 [Y 2 ]. (54) 

This system has two conservation laws: 

Stot = [So] + [Si] + [S 2 ] + [Yi] + [Y 2 ], E tot = [E] + [Yi] + [Y 2 \. 



From the steady-state equations derived from ( 53 ) and ( 54 ) we have that 

[Yi] = ^[Sm, [Y 2 ] = ^[S 2 ][E]. 

These expressions are increasing in [Si], [S 2 ], respectively. Using the conservation law for 
E to t and the two expressions for [Yi], [Y 2 ], we obtain 

^ k 7 k 9 E tot 



k 7 k 9 + k 6 k 9 [Si] + k 7 k 8 [S 2 ] ' 
which is positive and decreasing in both [Si] and [S 2 \. 
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Using the steady-state equation corresponding to (49) we have 

ks[Si] 



[So] 



(k 1 + k 3 )[EY 



which after substitution of [E] with the expression for [E] above, is increasing in [Si] and 

m 



We finally use the steady-state equation corresponding to (51)+(54) to obtain 



, r q i = k 2 k 7 k 9 E tot [S 1 ] k 3 k 5 

4[ 2j " k 7 k 9 + k 6 k 9 [Si] + k 7 k 8 [S 2 ] + h + k 3 [ lJ ' 

Fix a value of [Si]. The left-hand side of this equality is the line through the origin with 
slope ki in [S 2 ] ■ The right-hand side is a positive decreasing function of [S 2 ] defined for all 
positive values of [S2] ■ The function takes a positive value for [S2] = 0. It follows that the 
expressions on the two sides of the equality intersect in exactly one point for each fixed 
[Si]. This is the steady-state value of [S 2 ] corresponding to a given [Si]. Additionally, the 
left curve is independent of [Si] while the right curve increases in [Si]. As a consequence, 
[S2] increases as a function of [Si]. 

Using the remaining total amount Stot we have that 

Stot = [S ] + [Si] + [S 2 ] + [Yi] + [Y 2 ], 

where the right-hand side is expressed as an increasing positive function ip of [Si], such 
that ip(Q) = and such that it tends to infinity as [Si] does. Hence, for every value Stot > 
there exists a unique value [Si] > satisfying Stot = f(Si). Using this value of [Si], all the 
other steady states concentrations are positive and can be found using the relations above. 

We conclude that this model has exactly one positive steady state for all choices of rate 
constants and total amounts. 



4 Graphical illustration of the changes in number of steady 
states 

We illustrate here how the differences in number of possible positive steady states in the 
systems in Figure 2 arise from the differences in the conservation law shape. We choose the 
two-site phosphorylation representation as in the previous section with dephosphorylation 
reactions included. Since we would like to illustrate it in two dimensions, we cannot have 
two conservation laws. Therefore, we add the reactions 

-> S 2 , S 2 

to the system. The steady state behavior is not affected by this addition, but these reactions 
imply that there is only one conservation law, namely that of the kinase Etot- 

We use mass-action to model the core model and the extension model in column 1 of 
Figure 2 (with the modifications stated above). At steady state, the three models have the 
relation 

[So] = 7^2 ( 55 ) 
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for some constant Ai > 0. This relation is obtained from the steady-state equations and 
all models must have the same relation (Fact [3]). The other concentrations at steady state 
are derived from [E] and [So]. 



Therefore, the specific steady-state values are found by intersecting the graph of (55) 
with the curve given by the conservation law Etot- In the core model the curve is [E] = Etot, 
and hence it is a vertical line in the ([E], [So])-plane. For the extension models in the class 
in the first column of Figure XX, the curve is given by 



E tot = [E} + X 2 [E][S 



[So] 



E, 



tot 



[E] 



HE] 



Here A2 > is a constant that depends on the rate constants of the extension model. The 
steady states are positive if [E] < Etot- Note that the term that is added to the conservation 
law is different in the two models, because they belong to different steady-state classes. 
The form of this equation characterizes the class. 

We plot the steady-state curve for Ai = 2.5 (red) together with the conservation law 
equation with total amount E to t = 5 for two systems: green (core, corresponding to A2 = 0) 
and blue (extension model with A2 = 1). 




Steady-state curve 

1 steady state 

2 steady states 



For the core model and a particular choice of E to t the steady states are obtained by 
intersecting the steady-state curve (red) with the vertical line Etot = [E] (purple). For all 
values of E to t, the intersection has one point, hence there is one positive steady state. For 
the extension model, the steady states are found as the intersection of the steady-state 
curve with the expression for the conservation law Etot = [E] + A2 [-E 1 ] [So] (blue curve). For 
some values of A2, the two curves intersect in two points. 
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